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ABSTRACT 

Molecular hydrogen transitions in quasar spectra can be used to constrain variation in the 
proton-to-electron mass ratio, /i = m p /m e , at high redshifts (z > 2). We present here an anal- 
ysis of a new spectrum of the quasar Q0528— 250, obtained on VLT/UVES (the Ultraviolet 
and Visual Echelle Spectrograph, on the Very Large Telescope), and analyse the well-known 
H2 absorber at z — 2.811 in this spectrum. For the first time we detect HD (deuterated molec- 
ular hydrogen) in this system with a column density of log 10 (iV/cm~ 2 ) = 13.27 ± 0.07; 
HD is sensitive to variation in ft, and so we include it in our analysis. Using 76 H2 and 7 
HD transitions we constrain variation in /1 from the current laboratory value to be A/i//i = 
(0.3 ± 3.2 stat ± 1.9 sys ) x 10" 6 , which is consistent with no cosmological variation in /.i, as 
well as with previous results from other H2/HD absorbers. The main sources of systematic 
uncertainty relate to accurate wavelength calibration of the spectra and the re-dispersion of 
multiple telescope exposures onto the one pixel grid. 
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1 INTRODUCTION 



1.1 Molecular hydrogen 



The Standard Model of particle physics contains a number of "fun- 
damental constants", the values of which are predicted by the the- 
ory, and which must be measured by experiment. Two dimension- 
less constants of particular importance are the fine-structure con- 
stant, a = e 2 / (Aittahc), and the proton-to-electron mass ratio, 
H = m p /m e . The former determines the strength of electromag- 
netism, whereas the latter is a measure of the relative strengths of 
the strong and electroweak scales. Detection of spatial or temporal 
variation in a fundamental constant would directly demonstrate the 
incompleteness of the Standard Model, and might guide efforts to 
develop post-Standard Model theories. 

Quasar absorption lines can be used to constrain variations in 
a and n across most of the observable universe. In particular, the 
wavelengths of certain transitions are sensitive to changes in the 
value of q or yt; if a or \x were different, the relative wavelengths 
of the spectral lines would differ from those observed in the labo- 
ratory, even after correction for the redshift of the absorbing cloud. 
By comparing measurements of absorption transitions in quasar ab- 
sorbers to precise laboratory measurements, one is able to constrain 
changes in a or \x. 



I Thompson] {1975} noted that the Lyman and Werner transitions 
of molecular hydrogen (with Ao < 1150A) can be used to con- 
strain variations in /i at z > 2, where these UV transitions are 
redshifted into the optical spectrum and therefore can be observed 
with ground based telescopes. In particular, one can measure the 
quantity A/i//i = (fi z — fio)//-io where \i z is the measurement of 
/x at a redshift z and no is the laboratory value of fi. 

If (j, were to vary, the wavelengths of the H2 transitions would 
vary by different amounts. The shift can be quantified as 



A i = A°(l + *) (l + Ki^\ , 



(1) 



where A^ is the rest wavelength of the transition, Xi is the observed 
wavelength and Ki is a 'sensitivity coefficient' determining the 
magnitude and sign of the effect ( Varshalovich & Levshakov 1993}. 
That is, 



Ki 



d In A,: 
dln^i 



(2) 
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z is the absorption redshift of the cloud (equivalently, the redshift of 
a transition with no sensitivity to a change in n). The determination 
of z is not degenerate with determination of A/i//i provided that 
two or more transitions with different Ki values are used. 

It is worth noting that, if we define a nominal observation red- 
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shift for a particular transition as Zj, then 



Avj 
c 



1 + 2 



(3) 



That is, the velocity shift of an individual line is linearly propor- 
tional to Aft/fi. is known as the 'reduced redshift'. The redshift 
z is that of a transition with Ki = 0. In practice, 2 must be mea- 
sured simultaneously with Q from a set of transitions (e.g. through 
regression methods) because it is unknown a priori. 

Recent indications have emerged for a significantly positive 
A/i/jj, from the study of two H2 absorbers at z — 2.595 and z — 
3.025 toward Q0405-443 and Q0347-383 respectively. |Ivanchik| 
|et al . I ( [2005 | l analysed high resolving power (R ~ 50, 000) 
VLT/UVES spectra of these absorbers, and reported a combined 
value of Aft/ft for the two absorbers of (30.5 ± 7.5) x 10~ 6 
and (16.5 ± 7.4) x 10 -6 depending on whether H2 laboratory 
wavelengths were u sed from|Abgrall et ai| ( fl993| > or |Philip et al.| 
(2004) respectively. Reinho ld et al.| ([2606) analysed the same two 
absorbers using improved laboratory wavelength values, and Ki 
values derived from a more accurate calculation. Using a table of 
observed line wavelengths from Iva nchik et al.| l |2005} , they de- 
rived Afi/ft values of (27.8 ± 8.8) x 10 _e for the absorber to- spectrum of iFoltz et al. 
ward Q0405-443, (20.6 ± 7.9) x 10~ 6 for the absorber toward 0.005. jVarshalovich & 
Q0347-383, and a combined value of (24.5 ± 5.9) x 10" 6 . 

|King et al.| ( |2008| > analysed the same raw spectra of 
Q0405— 443 and Q0347— 383 using a more comprehensive fitting 
method in which all transitions are fitted simultaneously with the 
addition of Aft/ ft as a free parameter in the fit. They used an im- 
proved flux extraction and, importantly, an improved wavelength 
calibration procedure (Mu rphy et al.| |2007) which should signifi- 
cantly reduce systematic errors. They found Aft/ ft = (8.2±7.4) x 
10" 6 and (10.1 ± 6.2) x 10" 6 for the absorbers in Q0347-383 
and Q0405— 443 respectively. They also analysed the z = 2.811 
absorber toward Q0528— 250, providing a constraint of Aft/ ft = 
(— 1.4±3.9) x 10~ 6 . A weighted mean of these three values yielded 
Aft/ ft = (2.6 ±3.0) x 10" 6 . 

|Wendt & Reimersj ( [2008) and |Thompson et al.| ((2009> also 
analysed the absorbers toward Q0405-443 and Q0347-383 us- 
ing the same UVES spectra, and derived results statistically con- 
sistent with Aft/ ft — 0. Wendt & Molaro (2011 1 re-analysed the 
existing Q0347— 383 exposures combined with exposures under 
program ID 68. B-01 15(A) that had not been previously analysed, 
again deriving a result statistically consistent with zero. Given the 
result s of [King et al.|(|2008[>,|Wendt & Refiners] ( |2008| >, |Thompson| 
|et al.| (2009| > and |Wendt & Molaro| ( |2011fr it does not appear that 
there is significant evidence for Aft/ ft 7^ from the absorbers in 
Q0405-443 and Q0347-383. 

|Malec et al.| ( |2010| > reported a strong constraint from the 
z = 2.059 absorber toward J2123 — 0050 from spectra obtained 
using Keck/HIRES. A combination of good seeing conditions and 
the brightness of the quasar (r-band magnitude « 16.5 mag) meant 
that they were able to obtain a J? ~ 110,000 spectrum of suf- 
ficiently high SNR in a single night. They found that A ft /ft — 
(5.6 ± 5.5 stat ± 2.9 sys ) x 10 -6 from 86 H 2 transitions. Re- 
cently the same 12123—0050 absorber system was observed with 
VLT/UVES at a somewhat lower resolution (R ~ 53, 000) but 
with a higher signal-to-noise ratio (van Weerdenburg et al.|20lT| >. 
The spectrum covered 96 H2 and HD lines and yields a result of 
Aft/ ft = (8.5 ± 3.6 s tat ± 2.2 sys ) x 10 -6 , in very good agreement 
with the results from Keck. 

From consideration of the results above, it does not appear at 



present that there is any strong evidence for Aft/ ft 7^ at z > 1 
from measurements of H2 absorbers. 

Unfortunately the number of quasars known to contain H2 ab- 
sorption features in their spectra is low (< 15), which makes ob- 
taining a statistical sample of Aft/ ft measures from the H2 method 
difficult. To obtain a reliable measurement of Aft/ ft from H2, one 
needs a H2 column density that is high enough to produce signif- 
icant absorption, but preferably not saturated (which makes deter- 
mination of line centroids difficult). Similarly, transitions of signif- 
icantly different Ki must be used to achieve good precision on the 
measurement of Aft/ ft. Unfortunately, the only quasar absorbers to 
date which have been able to yield measurements of Aft/ ft at the 
< 10" 5 level are those listed above, which explains the significant 
number of analyses on these few objects. 



1.2 Historical Q0528 250 constraints 

Since the absorber in Q0528— 250 is the focus of this work, we note 
here previous constraints on variation in ft from this object. 

|Varshalovich & Levshakov|(|1993)) analysed the R ~ 4, 000 



1988f of Q0528-250 to obtain \Afi/fi 



'otekhin (1995} reanalysed the same spec- 



trum to conclude that \Aft/ft\ < 0.002. |Pote khin et al. ( 1998} used 
new observations of the same system at higher resolving power 
(R ~ 14, 000) to obtain Aft/ ft = (-10 ± 8) X 10~ 5 . 

|Cowie & Songaila| ( Tlt931 > used a R = 36, 000 Keck obser- 
vation of Q0528-250 to produce Aft/ ft G [-7, 5.5] x 10~ 4 (95 
percent confidence limits). 

Ubachs & Reinhold ( 2004 ) used the improved laboratory H2 
wavelengths o f|Philip et al.| (2004 ) to re-investigate the reported line 
positions for transitions in the spectra of Q0528— 250, Q0347— 383 
and Q1232+082. For the combined data, they found that Aft/ ft = 
(-0.5 ± 1.8) x 10" 5 . Omitting the Q0528-250 data, which was 
of poorer quality, they obtained Aft/ ft = (1.9 ± 1.5) x 10 -5 . An 
explicit value for Q0528— 250 was not given. 



1.3 Other quasar constraints on Aft/ ft 

The inversion transitions of ammonia are strongly sensitive to vari- 
ation in ft, with Ki ~ 4.2 ( |Flambaum & Kozlov||2007| >. |Mur-| 
|phy et ah] l |2008a| > and Henke et al. (20091 compared the inver- 



sion transitions of NH3 with rotational molecules to obtain strin- 
gent constraints on Aft/ ft at z < 1. Murphy et al. {2008a 1 used 
B0218+357 to obta in Aft/ ft = (0.74± 0.47 sta t ±0.76 5ys ) x 10" 6 
at z = 0.69, whilst |HenkeleTaTl ( |2009l used PKS 1830-2 11 to ob- 



x 10"° at z 



0.89. 



Kanekar 



tain A[i/ ft = (0.08 ± 0.47 ays 

(2011) analysed the z = 0.69 absorber toward B0218+357 to ob- 
tain Aft/ ft = (-3.5 ± 1.2) x 10~ 7 . 

Constraints derived from the analysis of H2/HD absorbers are 
not statistically competitive with the ammonia results at present. 
However, the ammonia results are from lower redshifts (z < 1), 
whilst the H2/HD results all probe z > 2, and therefore probe 
much larger distances and lookback times. The two methods are 

thus complementary. 

Recently it has been suggested that HaO + ions ( jKozlov 
|Levshakov|20 1 1 ) and CH3OH molecules ( |Jansen et al.|201i7 " 
even more sensitive probes to detect or constrain ft variation; these 
transitions exhibit large sensitivity coefficients (Ki > 10). 
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1.4 Motivation for re-analysis of Q0528 250 

The constraint on A/i//i derived from the Q0528— 250 absorber 
presented in |King et aL] (|2008 1 is very precise. However, there 
are good reasons to revisit this absorber. Some of the exposures 
which contributed to the spectrum used in that analysis are not 
well-calibrated. By well-calibrated, we mean that the exposures do 
not have ThAr (thorium-argon) calibration spectra taken immedi- 
ately afterwards. The design of VLT/UVES is such that the posi- 
tion of the spectrograph grating is reset between different quasar 
exposures. Although the specification is such that the placement of 
the grating should be good to within 0.1 pixels (D'Odorico et al. 
|2000| >, the use of spectra for which the ThAr spectra were taken 
after grating resets necessarily introduces wavelength calibration 
uncertainties into the spectrum. 

A more subtle concern for the earlier observations (under pro- 
gram IDs P66.A-0594, P68.A-0600 and P68.A-0106) is that they 
suffer from the fact that the slit width used for the observations was 
often significantly larger than the prevailing seeing conditions. For 
those observations, a 1 arcsecond slit was used in both the blue and 
red arms of UVES. The average ratio of the slit width to seeing, 
where seeing is quantified by the output of the DIMM (differen- 
tial image motion monitor) and the seeing values were weighted 
by the duration of the exposure, was 0.78. In one exposure of 1.6 
hours, the ratio was as low as 0.48. In the case where the seeing 
is significantly smaller than the slit width, the instrumental profile 
will be non-Gaussian, which complicates the analysis. Although in 
principle one can use a numerically-provided instrumental profile 
to convolve the Voigt profile model with, assuming that the instru- 
mental profile is Gaussian is significantly easier. 



2 DATA 

2.1 Description of the new Q0528 250 spectrum 

Note that this paper uses only new spectra of Q0528— 250 so as to 
measure A/i//i as independently as possible from that presented in 
|King et al.|p508) . 

Our new observations are summarized in Table Q] We ob- 
served Q0528-250 in late 2008/early 2009 under ESO program ID 
082.A-0087, with exposures totalling approximately 8.1 hours. All 
science exposures were followed immediately by ThAr calibration 
exposures without any intervening grating resets, thereby reducing 
drifts and shifts in the spectrograph in the intervening time. The av- 
erage slit- width-to-seeing ratio for these exposures is 1.03; that is, 
the slit width was more appropriately matched to the average see- 
ing conditions than for the exposures contributing to the spectrum 
analysed in King et al. ( 2008 ). 

We optimally extracted and wavelength calibrated the raw 
spectra using the ESO UVES Common Pipeline Language (CPL) 
software suite. While the CPL code redisperses the spectra onto a 
linear wavelength scale by default, we used only the original, 'un- 
redispersed' flux and error array for each echelle order in subse- 
quent reduction steps. A custom code, UVES_POPLEfQ was used 
to combine these extracted echelle orders from the 10 exposures 
into a single spectrum. The relevant heliocentric corrections and 
air-vacuum wavelength conversions were calculated and applied 
to the wavelength scale of each echelle order before redispersing 



1 See http://astronomy.swin. edu.au/~mmurphy/UVES_popler maintained 
by MIM. 



all orders onto a common, log-linear wavelength scale with disper- 
sion 2.0 km s -1 pixel -1 . The best-fitting relative flux scaling be- 
tween overlapping orders from all exposures was determined using 
an automatic \ 2 minimisation scheme. After applying these rel- 
ative scalings, the extracted quasar exposures were coadded with 
inverse-variance weighting and a cosmic ray rejection algorithm. 
Note that, as implied by equation [T| effects which shift the wave- 
length scale by a constant velocity at all wavelengths are generally 
unimportanj^] 

UVES_POPLER was also used to automatically fit a continuum 
to the final spectrum. This continuum was acceptable for our pur- 
poses only redwards of the Lyman-a forest. In the Lyman-a forest 
the continuum was manually re-fitted with low-order polynomials 
to obtain a nominal continuum against which the H2 absorption 
lines of interest could be reliably defined. Local constant or lin- 
ear modifications to this continuum were necessary in the vicinity 
of some H2 lines when conducting the detailed spectral fitting, as 
described below. 

Using the CPL pipeline we extracted the ThAr flux with the 
same spatial profile weights derived for the corresponding quasar 
exposure, thereby allowing the wavelength calibration polynomial 
most appropriate to the quasar exposure to be established. The CPL 
pipeline makes use of the ThAr line-list derived by Murp hy et al.| 
( |2007fr specifically for accurate calibration of UVES ThAr expo- 
sures. The calibration residuals had a root-mean-square (RMS) de- 
viation around the final polynomial of 70ms -1 in UVES's blue 
arm (where all the H2 lines fall) and 55 m s -1 over the majority of 
the two red arm chips, on average. 

From the widths of the extracted ThAr features we derived 
an average full-width-at-half-maximum (FWHM) resolution of 



^5.45 km s , as expected for a slit width of 0. 8 (D'Odorico et al. 



|2000j >. We initially used this FWHM value as the instrumental res- 
olution in our spectral fits. However, we expect - and later demon- 
strate - that a somewhat lower value of 5.15 km s -1 better repre- 
sents narrow, unresolved lines in the spectrum because the quasar 
illuminates the slit centrally rather than uniformly like the ThAr 
lamp. 



2.2 Laboratory wavelengths and values 

Until recently, the uncertainties in the laboratory wavelengths for 
the Lyman and Werner series H2 transitions were comparable 
to the typical uncertainties in the measured line positions in the 
quasar spectra. However, a significant amount of recent work has 
improved this situation by using laser-based spectroscopic tech- 
niques to obtain much more precise laboratory values. In partic- 
ular, |Philip et al.] pOOl) , |Ubachs & Reinhold| pOOl) and |Ivanov| 



et al. 



2008b 1 improve d the fractional w avelength accuracy (8X/X) 
whilst ! 



to ~ 5 x 10 
tional accuracy to 

-8 . 



Salumbides et al. 



1 2008 1 improved the frac- 



5 x 10 for m ost Lyman transitio ns and 



Bailly et al.|(|2009 



~ 10 for Werner transitions. See also 1 

The Ki coefficients can be calculated either through semi- 
empirical methods (Ubachs ct al. 2007) or from ab initio methods 
( |Meshkov et al.|2006 1. The Ki values derived from these two pro- 
cesses are in good agreement (to within 1 percent). 



2 This is not strictly true when many quasar exposures are combined, as is 
the case here. If different velocity shifts are applied to the different quasar 
exposures, and if the relative weights of the exposures vary with wavelength 
when forming the final, combined spectrum, then small relative velocity 
shifts will be measured between transitions at different wavelengths. 



4 J. A. King et al. 



Table 1. Journal of VLT/UVES observations of quasar Q0258— 250. Each row specifies the starting time of a single quasar exposure (column 1), all of which 
were taken with the DIC1 dichroic splitter, giving a central wavelength (A c ) in the blue ("B") and red ("R") arms as shown in column 2. The third column 
gives the exposure times for the quasar and subsequent ThAr wavelength calibration exposures in each arm. Column 4 gives the slit width used for each arm. 
These should be compared with the average FWHM seeing, as reported by the on-site Differential Image Motion Monitor (DIMM). The average difference 
between the spectrograph air temperatures at the times of the quasar and ThAr exposures, AT = Tqso — Tr-hAr* is given in the sixth column. The final 
column shows the atmospheric pressure difference between the quasar and ThAr exposures, AP = Pqso — ^"ThAr- On-chip binning by a factor of 2 in both 
spatial and spectral directions was used for all exposures. 



UT date & time 


Ac 




T ex p [S] 


Slit-width 


DIMM seeing 


AT 


AP 


[yyyy-mm-dd hh:mm] 


[nm] 


QSO 


ThAr 


[arc seconds] 


[arcseconds] 


[K] 


[hPa] 




B 


R 




B R 


B 


R 




B 


R 




2008-11-23 03:58 


390 


564 


2900.0 


1.6 1.8 


0.8 


0.7 


0.75 


0.2 


0.0 


0.2 


2008-11-23 04:52 


390 


580 


2900.0 


1.6 1.8 


0.8 


0.7 


0.79 


0.0 


0.0 


0.2 


2008-11-25 07:35 


390 


580 


2900.0 


1.6 1.8 


0.8 


0.7 


0.43 


0.0 


0.0 


0.0 


2008-12-23 05:03 


390 


580 


2900.0 


1.6 1.8 


0.8 


0.7 


0.75 


0.0 


0.0 


0.0 


2008-12-23 05:54 


390 


580 


2900.0 


1.6 1.8 


0.8 


0.7 


0.74 


0.0 


0.0 


0.1 


2008-12-23 06:46 


390 


580 


2900.0 


1.6 1.8 


0.8 


0.7 


0.99 


0.0 


0.0 


0.2 


2008-12-23 07:38 


390 


580 


2900.0 


1.6 1.8 


0.8 


0.7 


1.07 


0.1 


0.0 


0.1 


2009-01-25 02:40 


390 


580 


2900.0 


1.6 1.8 


0.8 


0.7 


0.72 


0.0 


0.1 


0.0 


2009-01-26 04:02 


390 


564 


2900.0 


1.6 1.8 


0.8 


0.7 


1.19 


0.0 


0.2 


0.1 


2009-02-26 02:52 


390 


580 


2900.0 


1.6 1.8 


0.8 


0.7 


0.75 


0.0 


0.0 


0.0 



|Malec et al.| ( |2010| > (see their Table 1) have collated the labo- 
ratory wavelength data from the above sources and the Ki values 
from |Ubachs et a l. ( 2007) to produce a comprehensive list which 
may be used for an analysis of A/i//i. They have used oscillator 
strengths from Abgrall et a l7| ([ 1 994) and calculated damping coef- 
ficients from |Abgrall et al."| i 2000 1. We therefore have used the lab- 
oratory wavelength data, Ki values, oscillator strengths and damp- 
ing coefficients tabulated in Malec et al. (2010) for our analysis of 
the H 2 absorber in Q0528-250. 



2.2.1 HD 

We have detected HD (deuterated molecular hydrogen) at z 
2.811 in the Q0528— 250 spectru m with a column density of 

log 10 (7V/cm- 2 ) 



Malec et al. 



(2010) have col- 



= 13.27 ±0.08 
lated oscillator strengths, laboratory wavelength values and Ki val- 
ues for HD; wavelength values are from Hollenstein et al. (2006) 
and |Ivanov et al.|(2008a| >, Ki values are from |Ivanov et al.| |2008a 
|2010[) and oscillator strengths were calculated by Malec et al 



2010) from Einstein A coefficients given in Abgrall & Roueff 



2006 1. We describe our treatment of HD further in section|3.5 



3 METHODS AND METHODOLOGY 

To estimate Afif/j, in Q0528-250, we model the H 2 /HD tran- 
sitions in the absorber with Voigt profiles using the non-linear 
least-squares \ 2 minimisation program VPFIT^] VPFIT was writ- 
ten specifically to fit Voigt profiles to quasar absorption spectra. 
Each Voigt profile is described by three parameters: the redshift of 
the transition, z, the column density, TV, and the Doppler width, 
b. Knowledge is also required of the rest wavelength, oscillator 
strength (/) and damping constant (IT 1 ) for each transition. The total 
model is the sum of the optical depths of a series of Voigt profiles 
convolved with a model for the instrumental profile (assumed to be 
Gaussian). 

If each transition in the absorber was well-represented by a 
single, unblended Voigt profile then determining Afi/fi would be 



http://www.ast.cam.ac.uk/~rfc/vpfit.html 



straightforward: one would measure a redshift for each H 2 transi- 
tion, and then calculate a set of Q values (equation [3}. Equation [3] 
implies that a plot of against would have a slope of A^i//i. 

In practice, there are two factors which complicate the anal- 
ysis. Firstly, the H 2 transitions all have rest wavelengths less than 
that of Lyman-a (1215.7A), and therefore are found only in the 
Lyman-Q forest, the dense series of broad H I transitions found 
bluewards of the quasar's Lyman-a emission line. This means that 
the H 2 transitions are seen against a background of relatively broad 
Lyman-a lines found at random redshifts and with varying opti- 
cal depths. To accurately determine the line centroids of the H 2 
transitions, one must model the Lyman-Q forest transitions simul- 
taneously with the H 2 transitions. Similarly, unless one allows the 
parameters of the Lyman-a transitions to be determined simultane- 
ously with those of the H 2 transitions, one will necessarily under- 
estimate the errors in the line centroids for the H 2 transitions and 
therefore also the error on A/j,/fi. 

Secondly, the Q0528— 250 absorber is not well-represented by 
a single Voigt profile. Instead, it shows 'velocity structure' - sev- 
eral clouds closely separated in velocity space but having different 
optical depths and Doppler widths. Provided that the SNR of the 
spectrum is sufficiently high, one can decompose the observed ab- 
sorption profile into different 'velocity components' (one for each 
absorption cloud). A reasonable model of the absorber is necessary 
if one is to obtain an accurate estimate of Afi/^i. 

In order to estimate A/j,/ n for the Q0528— 250 absorber it 
is therefore necessary to: i) construct a Voigt profile model for 
the Q0528— 250 H 2 absorber, determining the required number of 
velocity components using objective criteria, and; ii) construct a 
model for the Lyman-a forest in the vicinity of each of the H 2 
transitions. These two steps constitute almost all of the effort in 
accurately measuring Afi/ n from an existing spectrum of the ab- 
sorber (significant work being required to obtain a spectrum). In 
section [3~T1 we discuss the statistical criteria used to determine the 
preferred model. In section [3~2| we discuss how we build up a model 
of the Lyman-a forest in the vicinity of the H 2 transitions. In sec- 
tion |3.2.5] we describe how we determine the best-fitting number of 
H 2 velocity components. 
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3.1 Model selection 



Webb et al.| p999) |200T| |20TT) , |Murphy eTal] pOT) [20031 P0041 
2008b I and King et al. 2011) discussed the importance of accu- 
rately constructing a Voigt profile fit to an absorber in the context 
of potential variation of a. We apply a similar methodology here in 
attempting to create the simplest Voigt profile model which best ex- 
plains the absorption spectrum. King et al. (2008 ) and Malec et al. 
( |2010} discussed the importance of using an appropriate number of 
velocity components when fitting H2 absorbers to estimate Afi/fi. 

When comparing two models with the same number of free 
parameters, under the maximum likelihood method one should pre- 
fer whichever model has the lowest \ 2 - However, when comparing 
models with different numbers of free parameters, one must con- 
sider whether the reduction in \ 2 as a result of the addition of more 
free parameters is sufficiently large to justify the extra free parame- 
ters; a so-called 'information criterion' should be used to select the 
most appropriate model. 

One statistic which can be used to compare models with differ- 
ent numbers o f free parameters is the Akaike Information Criterion 
( Akai ke|"l974 1, defined as AIC = \ 2 + 2p, where p is the number 
of free parameters. As the AIC is only correct in the limit of large 
n/p (where n is the number of spectral points included in the fit), 
we use the AIC corrected for finite sample sizes, 
defined as 



AICC = x 2 + 2p+ 2 ^^ 
n — p — 1 



(4) 



The number of degrees of freedom, v, is therefore given by v — 
n — p. If several competing models are being considered, one 
chooses the model which has the lowest AICC. The actual value 
of the AICC is not important; only relative differences matter. The 
suggested interpretation scale for the AICC is the Jeffreys' scale 
(Jeffreys 196T) (seejLiddle 2007 for further discussion), where 
AAICC = 5 and AAICC = 10 are considered strong and very 
strong evidence against the weaker model respectively. 

When constructing the Voigt profile model, our goal is to find 
a physically plausible model for the spectral data which minimises 
the AICC. 



3.2 Constructing the Voigt profile model 

We built up our model of the molecular hydrogen transitions and 
the surrounding forest through an iterative process which we de- 
scribe here. With knowledge of the redshift of the molecular hy- 
drogen absorbers, in each spectrum we searched for molecular hy- 
drogen transitions which we considered to be potentially usable. 
We consider potentially usable transitions to be those for which the 
molecular hydrogen transition can be visually distinguished from 
its surrounds. This necessarily precludes the use of H2 transitions 
in regions of near zero flux, but in any event these transitions would 
contribute no meaningful constraint on Afi/fi. 



3.2. 1 Model construction process and physical assumptions 



We did not build up our region fits from those used in Ki ng et al.| 
(2008), but rather started afresh in order to provide as independent 
a measurement of Afi/ p as possible. From a list of potentially us- 
able transitions, we then selected a spectral region around the H2 
transition, where the region should be large enough to include any 
absorption feature which might overlap with the H2 transition. In 



general, we attempted to ensure that the fitting region was suffi- 
ciently large so as to return to the local continuum, although this 
was not always possible. In each of the fitting regions, we modelled 
the molecular hydrogen transition and then modelled all surround- 
ing features as H I. To do this, we added and removed H I compo- 
nents to attempt to achieve a statistically satisfactory model, using 
the criterion set out in section [3~T] The line parameters (N,b,z) for 
different H I transitions are independent of each other. Note that al- 
though most transitions observed in the forest are indeed due to H I, 
there are also metal transitions from other absorbers along the line 
of sight (including galactic and atmospheric lines). The identifica- 
tion of the origin of these transitions is not necessary if they do not 
overlap with the H2 transitions; we simply modelled them as H I in 
order to have a physical model for them. We describe the treatment 
of metal lines which overlap with H2 lines below. For all transitions 
assumed to be H I (which we refer to hereafter as just H I transi- 
tions), we use only the A1215.7 transition rather than the whole 
Lyman series, to prevent line misidentification spuriously impact- 
ing regions blueward of that transition. Where Lyman-/? transitions 
exist in the blue region of the spectrum, we simply modelled them 
with additional Lyman-a components. 

We then combined models from the regions fitted individu- 
ally into a model where the regions are fitted simultaneously. As 
the line parameters for the individual H2 transitions were indepen- 
dent when the regions were fitted independently, at this stage we 
imposed physical restrictions on the transitions by tying certain pa- 
rameters together to reflect the fact that the transitions are physi- 
cally connected. The ground states of the H2 and HD molecules 
have several rotational sub-states described by the quantum num- 
ber J. We refer to these sub-states as being different J-levels. For 
H2 and HD each J-level has a different ground state population. 
This means that the column density, N, should be the same for each 
velocity component arising from a particular J-level. The same ap- 
plies for the Doppler width, b, and the redshift, 2. We impose the 
assumption that the velocity structure is the same in all the J-levels, 
that is that a velocity component has the same value of z in all J- 
levels. We describe later the impact of this assumption. 

|King et aT] ( |2008| > did not actually impose the requirement 
that each velocity component arising from the same J-level has 
the same N; instead they fitted the column densities for each tran- 
sition as free parameters. We apply a similar approach. We must, 
however, ensure that a physical consistency is maintained, in that 
the ratios of the line strengths between different velocity compo- 
nents should be the same for transitions arising from the same J- 
level. We therefore imposed the requirement that the ratio of the 
column densities between the different components was the same 
for transitions arising from the same J-level. In this way, the total 
column density (effectively, oscillator strength) for each transition 
was a free parameter, but the ratios of the individual column den- 
sities of different velocity components within each transition were 
constrained. 

To model the forest, we thus iteratively refined the fit by alter- 
nately allowing VPFIT to minimise \ 2 f° r a particular model, then 
attempted to improve that model through the addition and deletion 
of H I components to obtain a robust model according to the criteria 
in sectionl3~71 



3.2.2 Blending with metal lines 

During the iterative process, it can become clear that a molecu- 
lar transition is blended with another line (presumed H I) when it 
was not thought to be from a fit to just that region. This is because 
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the information from the other molecular hydrogen transitions im- 
poses a strong constraint on the b parameters and redshifts of that 
transition, thus uncovering apparently hidden blends. These blends 
necessitate the addition of H I components that overlap with the H2 
or HD transition in question. With the addition of extra H I tran- 
sitions, an acceptable fit can generally be achieved. This demon- 
strates the utility of fitting all transitions simultaneously: otherwise- 
inconspicuous blends are generally revealed. 

In a few instances, the transitions which had to be included 
to achieve a statistically acceptable fit had extremely narrow b pa- 
rameters (6 < 5kms _1 ). In this case, it is likely that the blend is 
a metal line from an unknown absorber along the line of sight. As 
a result, we excluded such H2 or HD transitions from our fit. The 
reason for not accepting transitions affected by narrow-6 interlop- 
ers is that any inaccuracy in modelling the interloping transitions 
could lead to a significant bias in measuring the H2 line position - 
the narrow b parameter(s) of the interloping transition(s) means that 
the absorption they cause varies rapidly across the H2 line profile. 
Ultimately, the joint fit of all the molecular hydrogen transitions al- 
lows the detection and rejection of transitions which are likely to 
be contaminated by metal lines. 

Rejecting transitions which are suspected to be contaminated 
cannot bias Afi/fj, away from zero. Moreover, this should not bias 
Afi/fj. significantly. If the suspicion of contamination in particular 
lines was in fact due to A/i//i 7^ 0, we would expect to see this 
problem more frequently, and more obviously, for transitions with 
larger \Ki\. The number of transitions rejected was small, and did 
not appear to be correlated with \Ki\, and hence it is unlikely that 
we are biasing Afi/fj, towards zero. 



3.2.3 Continuum and zero level fitting 

Determination of the local quasar continuum is difficult in 
the Lyman-a forest due to absorption which can cover many 
Angstroms. Malec et al. (2010) discussed the fact that different in- 
vestigators may manually fit qualitatively similar but quantitatively 
different continuua to particular regions of a spectrum. For this rea- 
son, in any region where it appears that the local continuum has 
not been well determined, we allow for a freely floating continuum 
that is described by either a constant offset or a linear function, 
depending on the region in question. The use of a continuum fit 
which is at most linear helps to prevent strong degeneracies with 
the parameters of the Voigt profiles which might occur if a higher 
degree polynomial was used. The parameters of the continuum fit 
are determined simultaneously with all other parameters in the fit, 
and therefore any uncertainty in determining the continuum level 
naturally flows into the uncertainty on Afj,/fj,. 

Although weak night sky emission is subtracted as part of the 
initial flux extraction, the subtraction process does not appear to 
be optimal in the UVES CPL pipeline. In particular, in the base of 
saturated lines there appears to be a residual flux of about 2 percent 
of the local continuum. Therefore, in any region which includes 
saturated lines we allow for the zero level to be a free parameter in 
the fit, to be determined simultaneously with all other parameters. 



3.2.4 Over-fitting 

It is possible to add too many H I components to a particular region, 
leading to "over-fitting". Over-fitting is undesirable for several rea- 
sons. The primary reason is that it means that another, simpler 
model can explain the absorption spectrum better than the over- 



fitted model. Perhaps more importantly, it means that the perfor- 
mance of the optimisation algorithm can be substantially impaired. 
With significant over-fitting, convergence to the \ 2 minimum can 
be excessively slow. In extreme cases, convergence may not occur 
at all. Over-fitting can be detected through two means: 

i) The addition of components which increase the AICC sug- 
gests that the components are not supported by the data. If the 
AICC significantly decreases upon removal of the components, this 
suggests that the model was over-fitted. 

ii) Over-fitting causes the uncertainty estimates on the param- 
eters of the components in question to be excessively large (Gill| 
|et al.| [l986). In fact, this is often a good way to directly iden- 
tify components which are potentially unnecessary; the AICC re- 
lates to the model as a whole and therefore cannot suggest which 
components may be unnecessary. In particular, H I transitions with 
°"iog 1Q n > 1-0 or at/b > 1 are certainly suspicious. In regions 
with substantial over-fitting, errors can easily be substantially larger 
than this. The numerical cause of these large errors is strong rel- 
ative degeneracies between parameters. That is, \ 2 is almost flat 
in some direction in the parameter space relating to the offending 
transitions. It is this flatness in \ 2 which is the cause of poor con- 
vergence. Nevertheless, the presence of large errors on some com- 
ponents does not mean that they are unnecessary. In particular, the 
column densities for transitions which are saturated can be very 
poorly determined. This necessarily means that saturated H I tran- 
sitions will have large errors on the column density. 

Because of the impact of over-fitting on the convergence of 
VPFIT, we spent considerable effort trying to identify cases of over- 
fitting, and removed H I components as necessary to minimise the 
problem. 

3.2.5 Determination of the number of velocity components in the 
Q0528-250 absorber 

|King et ah] < |2008] > noted that at least 2 velocity components are 
plainly required to achieve an adequate fit to the Q0528— 250 ab- 
sorber. They argued on the basis of the AICC that 4 velocity com- 
ponents were required. We followed a similar procedure here, by 
comparing models with two, three and four velocity components. 
We choose whichever model has the lowest AICC (see section |4~2| 
below). 

3.3 Voigt profile fits 

Our final fits were obtained where we were not able to obtain any 
statistically appreciable improvement. We show one part of the 
H2/HD and Lyman-a Voigt profile fit in Fig.^ the full model may 
be found in the online-only version in Appendix [A] 

3.4 Determining A/i/fi 

There are two methods in the literature used for determining Afi/fj, 
for a particular H2 absorber. These can be described as the 'line- 
by-line fitting method' (LBLFM) and the 'comprehensive fitting 
method' (CFM). 

In the LBLFM, for each transition one calculates the £j values 
in equation [3] Equation [3] then implies that a plot of vs Ki will 
have a slope of Afi/fi. The slope of a linear fit to vs Ki can 
easily be obtained from standard \ 2 minimisation techniques. This 
method has the advantage of assisting visual detection of problem- 
atic transitions, where C,% deviates from the general trend shown by 
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Figure 1. Part of the H2/HD Voigt profile fit for the z = 2.811 absorber toward Q0528— 250. The vertical axis shows normalised flux. The model fitted to the 
spectrum is shown in green. Red tick marks indicate the position of H2/HD components, whilst blue tick marks indicate the position of blending transitions 
(presumed to be Lyman-o). Normalised residuals (i.e. [data - model]/error) are plotted above the spectrum between the orange bands, which represent ±lcr. 
Labels for the H2 transitions are plotted below the data. Other regions shown in Appendix[A](online only). 



8 J. A. King et al. 



other lines. A strong deviance of this nature implies that a problem 
exists in the fit to the transition which generated that Q. value; one 
can then return to the spectral fit to attempt to determine the prob- 
lem. This method was used by |Ivanchik et al.| l |2005| > and Reinhold 
[eTaLlpOOo] . 

Unfortunately, the LBLFM cannot be used when the absorber 
contains two or more closely spaced velocity components. The If 
velocity components overlap then the calculated errors in the line 
centroids for those velocity components will be correlated. It fol- 
lows that, to use the LBLFM, one must first demonstrate that a 
single component is the best model for the H2/HD lines i.e. that 
a two component model is statistically worse than a one compo- 
nent model. The LBLFM can be generalised to fit closely spaced 
velocity components through the use of generalised least squares 
methods, which allow for correlations between data points. How- 
ever, to the best of our knowledge this has not been used in the 
literature to determine Afi/fi from a quasar spectrum. 

The CFM allows one to reliably measure A^i//i in the case 
where the absorber contains two or more closely spaced velocity 
components. In the CFM, one fits all H2 transitions simultaneously 
with a single redshift for each velocity component, and allows for 
A/i/jj, as a free parameter in the fit, such that the rest wavelengths 
are perturbed as \® — > A°[l + Ki(Afi/ fi)]. The value of Afi/n at 
the x 2 minimum is the best estimate of Afi/n from the spectrum. 
In this way, the overlap of the velocity components is naturally ac- 
counted for within the fitting process, including the error estimate 
for A/1//1. Additionally, this method significantly reduces the num- 
ber of free parameters in the fit, which should improve the reliabil- 
ity of the method. This method was used by |King et alT]p008) and 
|Malec et al.|pOTo) . 

We have assumed that the errors in the Ki values (which are 
~ 1 percent) are negligible. This assumption can be clearly justified 
by examining typical reduced redshift plots e.g. Fig. 1 in ( |King| 
|et al.|20 08 ), where it is clear that the uncertainty in determining the 
velocity shifts of the line centroids is much more important than 
~ 1 percent errors in the K, coefficients. 

The absorber in Q0528— 250 contains several closely spaced 
molecular velocity components. We therefore use the CFM to de- 
rive our estimate of A/1//1. We note that our fits to each region were 
all constructed assuming A/i//i = 0; A/i//i was allowed to vary 
only once our model for the H2/HD and Lyman-Q lines was com- 
plete. If Afi/fj, is significantly different from zero, it is conceivable 
that this method could bias our measurement of A/1//1 towards zero 
by "fitting away" any /i variation through the choice of the Lyman- 
a forest model. However, because the forest lines generally have 
high b parameters and because our model cannot fit the spectra in 
a totally arbitrary fashion (our model consists of a series of Voigt 
profiles, which have rather less flexibility to fit the spectrum than 
polynomials, for example) any bias induced should be small. Note 
that, on account of the above methodology, we cannot bias Afi/fi 
away from zero. 

3.5 HD 

HD is sensitive to a change in fj,, and therefore here we include HD 
in our analysis of Afi/fi. Although HD should display a similar 
velocity structure to H2, the low optical depth of the HD transitions 
and the small number of transitions observed means that any such 
structure is unresolved. We therefore model the HD absorption with 
only a single velocity component. That is, the constraint on Afi/fi 
from HD is derived only by considering potential velocity shifts 
of the HD lines with respect to each other, and not with respect 



to H2 transitions. The small number of HD transitions used means 
that the statistical constraint on Afi/fi derived from HD is weak 
compared to that derived from the H2 transitions. However, future 
detected HD absorbers observed with greater optical depths and 
higher SNRs may provide statistically competitive constraints. 



3.6 Spectral extraction problems 

The CPL pipeline appears to incorrectly estimate the errors associ- 
ated with the flux data points in the base of saturated lines. In partic- 
ular, the dispersion of the flux data points about the local zero level 
is too large to be accounted for by the statistical error. The statistical 
errors on the flux points appear to be under-estimated by a factor of 
at least 2. Although it is difficult to determine precisely what hap- 
pens in regions of low, but non-zero flux, we believe that the errors 
there are also underestimated. The effect of this is to give falsely 
high precision on any quantity derived from these data points (in- 
cluding Afi/fjb), Additionally, one cannot fit plausible models to 
data involving regions of low or negligible flux. 

When we combine the individual exposures into a final spec- 
trum using UVES_POPLER, the program calculates a check on the 
consistency of the exposures contributing to each pixel. As the 
combined value of each flux pixel is given as a weighted mean 
of pixels from the contributing spectra, this consistency check is a 
value of xit for each pixel about the weighted mean. We attempt to 
correct for the under-estimation of uncertainties in regions of low or 
negligible flux by applying the following algorithm. For each pixel 
in the combined spectrum, we take a region of five pixels centred 
on that point. We then take the median of the %v values just de- 
scribed that are associated with those five points. We then multiply 
the error estimate for that spectral pixel by the square root of that 
median value (that is, en — > <n x \J median[xS])- 

However, we only increase the error array if median(xS) > 
1. This is because the spectral errors are generated from photon- 
counting, and therefore the statistical uncertaintiy should be a lower 
bound on the true uncertainty. Thus, it is difficult to justify decreas- 
ing the error array without particular evidence that the spectral un- 
certainty is systematically under-estimated. 



4 RESULTS 

In section |4~2| we describe our determination of the velocity struc- 
ture of the absorber, and the A/1//1 values (and associated statistical 
errors) which result. In section|43]we describe various consistency 
checks we make to ensure that our results are robust. In section [5] 
we detail our analysis of potential systematic errors. In section|6]we 
give our preferred result including the final estimate of systematic 
errors. 



4.1 Transitions used 

For our analysis, we have included 76 molecular hydrogen transi- 
tions. In Fig.|2]we give the relationship of Ki and J with Ao for the 
H2 and HD transitions used in our analysis. The seven HD transi- 
tions we have used are: L1R0, L2R0, L3R0, L4R0, L5R0, L6R0, 
L8R0. We show the spectral regions for the stronger HD transitions 
as a velocity plot in Fig. [3] clearly demonstrating the existence of 
HD. 
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Figure 2. Relationship of K t and J with Aq (rest wavelength) for the H2 
and HD transitions used in our analysis of the 2 = 2.811 absorber toward 
Q0528-250. Only H 2 transitions with J e [0, 4] are shown. All the HD 
transitions have J = 0. Upper panel: the sensitivity coefficients, Ki, for 
the transitions used in our analysis of the absorber (dark blue points) and 
not detected or not fitted (grey points). Lower panel: the distribution of 
transitions with wavelength according to their J-level. 



4.2 Velocity structure and Afi/fi results 

We examine the question of the velocity structure of the 
Q0528— 280 absorber under different scenarios. To do this, we con- 
sider models with two, three and four velocity components. In the 
3- and 4-component case, we applied a model in which correspond- 
ing velocity components in transitions from the same J-level had 
the same b parameter [b — F(J)], and also the scenario in which 
corresponding transitions had the same b parameter regardless of 
J [& / F(J)]. For the 2-component case, we only considered 
b = F(J). We give the AICC for these scenarios, and the result- 
ing values of Aft/fj,, in Table [2] We note that we were unable to 
obtain a stable fit for a 4-component model where b parameters 
for corresponding components were forced to be the same for all 
J-levels. In this model, the column density of one of the compo- 
nents in the J = 1 transitions was driven down to effectively zero. 
This component was the second strongest component in the J — 2 
and J = 3 transitions. Although we could have omitted this com- 
ponent, the substantial differences in relative strength between the 
different components in the different J-levels means that this model 
is very unlikely to be a good representation of the physical situa- 
tion, and therefore that the value of Afi/fj, derived is unlikely to be 
accurate. 

We note the sign change in A/i//i from the result in Kin g et af] 
1 2008| [A/i/V = (-1.4±3.9)xl(T 6 ], although in then = 3 case 
the result is only marginally different from zero. It is clear from the 
n — 3 results that a model where different J-levels have different 
b parameters is preferred very strongly over a model with the same 
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Figure 3. Different spectral regions for our spectrum of Q0528— 250 
demonstrating the presence of HD in the absorber at 2 = 2.811. The veloc- 
ity is shown about 2 = 2.81112, which is the fitted HD redshift. The black 
line is the spectral data and the green line is the model fitted to the data. 
The red tick marks at V = are fitted HD transitions. The red tick marks 
at v 7^ are fitted H2 components. The blue tick marks are fitted Lyman-a 
components. 



Table 2. Analysis of the velocity structure of Q0528— 250. n gives the num- 
ber of components. The second column defines whether the 6-parameter 
for different components is fixed or is different for transitions with differ- 
ent J, i.e. whether 6 is a function of J [b = F( J)]. The column labelled 
AICC gives the AICC for the fit. The column AAICC shows the difference 
of the AICC with respect to the best-fitting model. For the 4-component, 
b ^ F( J) model, the column density of one component of the J = 1 tran- 
sitions was driven to zero, and the component was rejected, where this com- 
ponent was strongly detected in other J-levels. This implies that the model 
is not physically realistic, and so we label it as unstable. For these fits, the 
instrumental resolution was assumed to have a FWHM of 5.45 kms -1 . 
Note that the AICC is numerically quite large because of the thousands of 
degrees of freedom due to thousands of spectral points being included in 
the fit. 
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b = F(J) 


AICC 
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A M / M (10- 6 ) 
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b = F(J) 
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b = F(J) 


11488.2 





1.115 


0.2 ±3.2 


3 


b # F( J) 
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1.141 


0.4 ±3.2 
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b = F(J) 
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3.4 ±3.7 
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Unstable fit 
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6-parameter for each J-level. This accords well with the results of 
|King et al.| < [2008) . 

We note that, for these fits, the 3-component model is pre- 
ferred to the 4-component model, at about the same statistical 
significance as the 4-component model was preferred to the 3- 
component model in King et al. (2008). There are several points 
to note here: 

i) The seeing conditions for the spectra used by |King et aL] 
( 2008 1 were quite variable, and may have induced a significantly 
non-Gaussian instrumental profile. The requirement for 4 compo- 
nents in King et al. (20081 may be a reflection of the non-Gaussian 
profile, rather than the absorber itself. Because the components are 
unresolved at these resolving powers, identification of the correct 
number of components is difficult. Nevertheless, what is important 
is that all the statistical structure in the measured profiles is ade- 
quately accounted for. This is directly reflected in the AICC. 

ii) Some of the exposures contributing to the spectrum used 
by |King et aT] l |2008l were poorly calibrated. It is conceivable that 
wavelength miscalibrations could cause a model with more com- 
plexity to be favoured. 

Hi) Although we have attempted to apply the same forest 
model in analysing the 3- and 4-component model (where in each 
case x 2 is obviously minimised with respect to all the parameters), 
the construction of the forest model itself depends on the choice 
of the H2 model. Strong H2 components will obviously have little 
effect on the forest model because they are clearly distinguished 
from the forest. However, the 4th component of the model is weak 
and unresolved visually. We created the 3-component fit by remov- 
ing the weakest component from the 4-component fit. In most re- 
gions, the resulting fit was reasonable, however in a small number 
of regions we found that we had to add weak forest components to 
account for the removal of the 4th H2 component. To achieve a like- 
with-like comparison, we included these extra forest components in 
the 4-component fit. This means that any test for the statistical sig- 
nificance of the number of components depends somewhat on the 
choice of forest model near the H2 lines. 

iv) We noted whilst we were iteratively refining the model in 
the 3- and 4-component cases that the 4-component model was pre- 
ferred for most of the refining process, with AAICC ~ 10 in its 
favour for much of the time. It was only in the last few rounds of 
refining the model that the 3-component model became preferred as 
a result of changes made to a small number of regions. Therefore, 
the choice of the correct number of components can be sensitive 
to decisions made about the forest model in a small number of re- 
gions. 

v) If the instrumental resolution used differs significantly from 
the true instrumental resolution then one might choose an incorrect 
number of components. 

We conclude from this that although the Jeffreys' scale sug- 
gests that there is very significant evidence for the 3-component 
model over the 4-component model, in light of the fact that this 
evidence is conditioned on the correct choice of forest model and 
instrumental resolution, the actual preference for the 3-component 
model over the 4-component is rather weak. These arguments ap- 
ply similarly to the preference for a 4-component model over a 3- 
component model in King et al. (2008). The issue of whether there 
are 3 or 4 components is simply very difficult to resolve given the 
actual SNR and resolution of the spectra available. 

To investigate the effect of varying the instrumental resolu- 
tion, we re-ran the 3- and 4-component, b = F( J) fits with an in- 
strumental resolution of 5.15 kms -1 . The results of this are given 
in Table [3] There are three points to note from varying the instru- 



Table 3. Effect of varying the instrumental resolution from 
FWHM=5.45kms _1 (table [5} to FWHM=5.15 km s _1 on estimates 
of Afi/fi for different models, n gives the number of components. The 
second column defines whether the fc-parameter for different components 
is fixed or is different for transitions with different J, i.e. whether b is a 
function of J [b = F(J)]. The column labelled AICC gives the AICC for 
the fit. The column AAICC shows the difference of the AICC with respect 
to the best-fitting model. 



n 


? 

6 = 


F(J) 


AICC 


AAICC 


A M / M (10- 6 ) 


3 


b = 


F(J) 


11470.4 





0.3 ±3.2 


4 


b = 


F(J) 


11500.1 


29.8 


4.6 ±3.8 



mental resolution: i) in both cases the AICC decreases, suggesting 
that the lower instrumental resolution is preferred; ii) the preference 
for the 3-component model over the 4-component model increases, 
and; Hi) A/i/ [i is relatively insensitive to the instrumental resolu- 
tion, as one would naively expect (particularly for the 3-component 
model). We give the resulting parameters from this fit in Table [4] 

In principle, one can account for the uncertainty in the choice 
of model by constructing a weighted mean of the values 
from each model, where the weights are given by the penalised like- 
lihood [i.e. Wi — exp(— AICC/2)]. Because the difference in the 
AICC between the 3-component and 4-component models is « 30, 
however, the suppression of the 4-component model is sufficiently 
large so that any contribution from this model can be neglected. 

Thus, on the basis of the statistical results in Table [3] we 
choose Afi/fi — (0.3 ± 3.2) x 10 -6 as our preferred statistical 
result for the analysis of z — 2.811 absorber toward Q0528— 250. 
This is obviously consistent with no change in p. 



4.3 Consistency checks 

We can relax some of the assumptions made in our analysis of this 
absorber to explore whether they have a meaningful impact on the 
result. In particular, we explore here whether the result we obtain 
is significantly affected by our assumptions about the different J- 
levels. We follow a similar procedure to that used by Mal ec et al.| 
l |2010) . We also verify that our optimisation algorithm is function- 
ing adequately. 



4. 3. 1 Different A^i//i from different J -levels 

Rather than allowing transitions from all J-levels to contribute to 
a single value of Afi/fi, within VPFIT we can calculate a value of 
A/j,/ n for each H2 J-level (and one for HD) separately. Strictly, 
the values of A/1//1 obtained are not independent because: i) they 
assume that each J-level has the same number of velocity compo- 
nents; and, ii) the redshifts are tied between corresponding velocity 
components in transitions arising from different J-levels. Never- 
theless, this is useful for quantifying the contribution that each J- 
level makes to the final result. Ubachs et al. (2007 ) noted that, on 
account of the para-ortho distribution of H2 the J = 1 state is 
significantly populated even at low temperatures. They suggested 
dividing the states into a J € [0, 1] set (cold states) and J ^ 2 set 
(warm states) to examine the impact of temperature. We examine 
both of these cases in Fig. [4] We see that there is no clear evidence 
for a difference of Afi/fi obtained using transitions arising from 
different J-levels. 
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Table 4. Best-fitting H2/HD model parameters in the 3 velocity component model by J-level. n gives the number of H2 transitions fitted for that J-level. 
The columns labeled "Component 1" through "Component 3" give the parameters for the 3 H2 components, whilst the HD column gives the parameters 
for the single fitted HD component. For the H2 components, iV re i specifies the relative strengths of the three components, defined as AT re x j = NifEjNj 
for the components in that J-level; because we fit the line intensities as free parameters for the H2 transitions it is difficult to derive an accurate estimate 
of the column density. For an approximate estimate of the total column density for each J-level, see Table [6] The HD components were fitted using the 
actual oscillator strengths, and so we give a direct estimate of the column density for the J = level of HD explicitly, b gives the b parameter for the 
component. Uncertainty estimates are derived from the covariance matrix at the purported optimisation solution. For certain poorly-determined b parameters, 
where cr^/b > 1 these uncertainty estimates will not yield accurate confidence intervals on the b-parameters due to the fact that b > 0. Our optimiser has a 
hard lower limit for b of 0.4 km s" 1 to prevent numerical problems with the convolution of the Voigt profile model with the instrumental profile model, which 
is why the b parameters for the J = 3 and J = 4 levels of component 3 are both 0.40 km s~ 1 . Note the trend for decreasing b with increasing J for the H2 
components, which is why the b = F(J) models are significantly preferred over the b ^ F[J) models. 



J-level 


n 


Parameter 


Component 1 


Component 2 


Component 3 


HD 






z 


2.8110056(9) 


2.8111223(7) 


2.8109346(11) 


2.8111200(33) 


J = 




b (kms- 1 ) 


5.36 ±0.36 


3.56 ±0.22 


1.09 ±0.94 


0.95 ±0.87 


3 




0.68 


0.25 


0.08 


log 10 (iV/cm- 2 ) = 13.27 ± 0.07 


J = 1 




b (kms- 1 ) 


4.64 ± 0.99 


3.57 ± 0.25 


2.44 ± 0.62 




8 




0.69 


0.22 


0.09 




J = 2 


24 


b (kms- 1 ) 


4.63 ±0.25 
0.18 


3.05 ±0.08 
0.07 


0.42 ±0.08 
0.75 




J = 3 


28 


b (kms- 1 ) 


3.78 ±0.19 
0.14 


3.08 ± 0.08 
0.07 


0.40 ± 0.05 
0.78 




J = 4 


13 


b (kms -1 ) 


3.58 ±0.26 
0.42 


1.88 ±0.15 
0.46 


0.40 ±0.22 
0.12 
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Figure 4. Left panel: A/1//1 for each H2 J-level and HD, assuming that 
each J-level has the same 3-component velocity structure. The J = level 
has a substantially larger uncertainty than the J £ [1,4] transitions because 
only 3 J = transitions are used in the fit. HD has been plotted in units 
of 10 — 5 to increase clarity for the H2 results. Right panel: A.p/ p for two 
groups of transitions, H2 J £ [0, 1] ("cold transitions") and J S [2, 4] 
("warm transitions"). The black data points are the results where the red- 
shifts for corresponding components in the cold and warm transitions are 
assumed to be the same, and the red points are where the redshifts are al- 
lowed to differ between the cold and warm transitions. 



4.3.2 Optimiser performance 

The determinisation of Ap/p requires minimising \ 2 for a model 
which contains hundreds of free parameters. Adequate perfor- 
mance of the optimiser is crucial in obtaining a reliable value of 
Ap/p and an associated error. Premature termination of the opti- 
miser will yield an erroneous best estimate for Ap/p. Typically 
this result will be biased towards the value of Ap/p that is used 
to initalise the optimiser. Code errors or other gross problems with 
the optimiser will also render the error associated with Ap/p un- 
reliable. 

In the context of potential variation of a, we have checked that 
the optimiser within VPFIT performs adequately for the case of a 



few tens of parameters using Markov Chain Monte Carlo (MCMC) 
methods l |King et "al7||2009| >. We found for three systems that the 
estimated value of Act I a and its associated la error provided by 
VPFIT are very similar to that determined from MCMC exploration 
of the likelihood function. Nevertheless, our model for the absorber 
in Q0528— 250 and the Lyman-a forest have more than an order of 
magnitude more free parameters, and so an explicit check on the 
performance of the optimiser is advantageous. 

A simple check on whether the optimiser is functioning ade- 
quately is to restart VPFIT with an initial value of Ap/p that is sig- 
nificantly different to both the original starting guess (Ap/p = 0) 
and the final result returned by VPFIT. To do this, we re-ran the 
optimisation with a starting value of Ap/p = 10 _J for the 3- 
component, b = F( J) fit with an instrumental resolution FWHM 
of 5.45 kms -1 . The final value of Ap/p under this circumstance 
was 0.210 x 10~ 6 , compared to 0.187 x 10 -6 when started from 
Ap/p — 0. These two numbers differ by ~ 0.007a, which is en- 
tirely negligible. This demonstrates both that the final result is in- 
sensitive to reasonable choices of the starting guess for A p/ p and 
that our optimiser is functioning adequately. 

We note that Malec et al. (2010) performed Monte Carlo sim- 
ulations on realisations of their J2123— 0050 spectrum and found 
consistency between the error on Ap/p and the distribution of 
Ap/p values. Similar simulations confirm the reliability of VPFIT 
in the context of Act /a l |Murphy|2002{ |Murphy et al.|2003||King| 
|et al.|2009| > 



5 SYSTEMATIC ERRORS 



Malec et al. (2010) discussed a number of potential systematic er- 
rors which could affect the measurement of Ap/p from the analy- 
sis of H 2 /HD absorbers. We investigate the same potential system- 
atic errors here. 

It is worth noting initially that a clear potential source of un- 
certainty arises from the quality of the wavelength calibration. Be- 
cause the Ki values are well-correlated with wavelength within the 
Lyman and Werner series (but much less so when considering both 
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the Lyman and Werner series simultaneously), any effect which 
systematically expands or compresses the wavelength scale will 
mimic a change in fj,, Any effect which only locally perturbs the 
wavelength calibration will have a random effect on Afif/j,, and, 
providing enough transitions are used, the impact of such an effect 
will average to zero. We discuss wavelength calibration in the next 
two sections. 

5.1 Known wavelength calibration errors due to 
uncertainties in the ThAr calibration 

The calibration of the ThAr wavelength scale is not perfect; each 
of the ThAr transitions displays a residual velocity offset about 
the best-fit polynomial solution. The RMS of the residuals is ~ 
70 ms -1 in the blue arm and ~ 55 ms -1 in the red arm. How- 
ever, these fluctuations are random, and therefore will average out 
if a large number of H2 transitions are used. Only systematic devi- 
ations from the true wavelength solution should appreciably affect 
the best estimate of Afi/fi. There are fewer good ThAr lines in the 
blue end of the spectrum than in the red end, and therefore larger 
deviations of the wavelength solution from the true solution are 
possible. 

Following the analysis of Murphy et al. (20071, the system- 
atic deviation in the blue end of the spectrum relative to the red 
end of the forest has an upper limit of ~ 20 m s -1 . The maximum 
Ki value used in the fit is 0.053, whilst the minimum is —0.009. 
This implies that the maximum possible systematic due to this ef- 
fect is given by 5(Afi/fi) = (Av/c)/AKi, which is 1.1 x 10 -6 . 
In reality, the effect is likely to be smaller than this as positive de- 
viations should tend to cancel somewhat with negative deviations. 
However, how to reduce the effect is unclear; it may not simply 
scale as 1/vJV (for TV transitions). Therefore, we retain this esti- 
mate as the maximum possible systematic effect due to this cause. 

5.2 Intra-order distortions of unknown origin 

The path that the quasar light takes through the telescope is similar 
but not identical to that from the ThAr calibration lamp - the ThAr 
light fills the slit nearly uniformly, whilst the quasar light does not, 
and the ThAr light does not pass through the telescope optics. Due 
to the different light paths, the wavelength scale of the quasar light 
may be different to that of the ThAr light; the differences between 
them may appear as an apparent distortion of the wavelength scale. 
Both long range and short-range distortions are possible. 

|Griest et al.| (2010) identified a pattern of distortion within 
echelle orders in Keck/HIRES spectra, such that the wavelength 
scale at the centre of echelle orders is distorted with respect to that 
at the echelle order edges. The peak-to-peak velocity distortion is 
~ 500 ms - at ~ 5600A. The distortion was identified by com- 
paring the calibration of a spectrum obtained using a ThAr expo- 
sure to that obtained using an I2 absorption cell. The iodine cell 
is placed in the quasar light path, and the characteristic absorption 
spectrum is imprinted on the quasar spectrum. The use of an iodine 
cell therefore obviates the concern about optical path differences 
when using a ThAr lamp. Unfortunately, an I2 cell is not useful for 
calibration of general quasar observations, because the iodine tran- 
sitions cover only a relatively narrow part of the optical range, and 
because of the loss of flux from the quasar as a result of the use of 
the cell. The observed distortion pattern appears to be dependent 
on wavelength, and the distortion may be larger at longer wave- 
lengths. The precise origin of the distortions is unknown, and sim- 
ilarly it is unknown to what extent the distortions remain constant 



in time, and how they depend on extrinsic factors (e.g. telescope 
orientation, temperature, pressure and accuracy of quasar centering 
in the spectrograph slit). Therefore, it is not possible at present to 
adequately remove these distortions of the wavelength scale from 
observations. 

Whit more et al.] (2010| ) identified a similar effect in 
VLT/UVES spectra, with a peak-to-peak velocity distortion of 
~ 200 ms -1 . The distortion appears to be much less consistent 
between echelle orders than that seen by |Griest et al.[ however. Fur- 
ther observations by one of us (MTM) have shown that the observed 
wavelength distortion is definitely not constant over long periods of 
time i.e. more than several nights, which makes removal of the dis- 
tortion extremely difficult. 

Similar to Malec et al. (2010), we attempted to estimate the 
magnitude of the error introduced into a determination of Afi/ /i 
as a result of the observed velocity distortions. To do this, we 
used a triangular-shaped distortion, where the wavelengths of pix- 
els at the centre of echelle orders were displaced by +200 ms -1 
with respect to those at the echelle order edges. The modification 
to the spectra was implemented within UVES_POPLER. The shift 
in Afj,/ n after modifying the spectrum was —0.3 x 10 -6 [using 
the 3-component, b = F(J) model]. Clearly this value is model- 
dependent - if the distortion has a different amplitude or form, then 
the impact on A/j,/fi may be different. However, this estimate of 
the systematic error is likely to be of the correct magnitude. We 
therefore adopt a Gaussian with a — 0.3 x 10 -6 as an estimate of 
the systematic effect due to distortions of this type. 

5.3 Velocity structure and spatial segregation 

It is possible that transitions arising from different J-levels might 
be spatially segregated (Jenkins & Peimbert 1997, Levshakov et al. 
2002). Assuming that all J-levels arise from the same redshift in 
this event could spuriously produce A/i/^i 7^ 0. Although our 
result for the Q0528— 250 absorber is statistically consistent with 
zero, it is of course possible that a non-zero A^i//i could be pushed 
towards zero by this sort of systematic effect. Similar to Malec et al. 
l |2010) , we relaxed our assumption that corresponding components 
in all J-levels arise from the same redshift. In particular, we divided 
the data set into "cold" transitions, J G [0, 1], and "warm" transi- 
tions, J £ [2, 4], as was done earlier, but only tie the redshifts of 
corresponding components between different J-levels within these 
two groups. If there is spatial segregation, this should be seen as 
a statistically significant difference between the redshifts of corre- 
sponding components between the two groups, and also as a sub- 
stantial shift in the values of A^i//i derived from the two groups 
compared to what was obtained earlier. 

In the right panel of Fig. [4] we directly compared Afi/fj, in the 
case where the velocity structure was allowed to vary between cold 
and warm components, and note that there is no appreciable shift. 
In Fig.[3]we show these considerations more directly by examining 
the differences in the redshifts of the three components, and also the 
explicit difference between A/1//1 in the two cases considered. We 
see that there is no statistically significant difference between the 
redshifts in any of the three components when the velocity structure 
is allowed to vary between the cold and warm components. Simi- 
larly, we see that the shift in A[i/ fi is < 0.1a when the velocity 
structure is allowed to vary between the cold and warm compo- 
nents. Thus, we conclude that there is no evidence for a systematic 
shift in A[i/ fi as a result of segregation of the cold and warm J- 
levels. 

Nevertheless, to quantify the possible error introduced by our 
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Figure 5. Left panel: Difference between the redshifts of components for 
the H2 fit for the 3-component model to the absorber when the velocity 
structure was allowed to vary between the "cold" and "warm" components. 
The difference for each component is defined as z( J £ [0, 1]) — z( J £ 
[2, 4]). Right panel: The difference in Ap/p, when the velocity structure 
was allowed to vary, defined as A^t/^t(structure allowed to vary) — 
A/i//^(structure not allowed to vary). The error estimate is calculated 
as the mean of the error estimates in the two cases considered; the errors do 
not differ appreciably between the two cases. From consideration of both 
of these plots there does not appear to be any significant evidence that our 
assumptions about the velocity structure model have affected our estimate 
of A(j,//i. 



assumptions regarding velocity structure, we examine the actual 
shifts in Ap/ p. The shift in Ap,/ p for the J £ [0, 1] levels from 
the original value [for the 3-component, b = F(J) model] is 
0.8 x 10" 6 , and for the J € [2, 4] levels is 0.2 x 10~ 6 . We thus 
take 0.8 x 10 -6 as an estimate of the potential error introduced 
into our analysis due to assumptions about the velocity structure in 
order to be conservative. 



5.4 Re-dispersion of spectra 

The spectrum used is the result of the co-addition of exposures 
taken at different times, each with its own ThAr calibration. Dur- 
ing the co-addition, the spectra are placed on a common wavelength 
grid. Because the spectra are re-binned, the choice of wavelength 
grid introduces correlations between neighbouring pixels. More 
importantly, the choice of the wavelength grid has the potential to 
affect A/i/jit by slightly distorting the absorption profile shapes. To 
investigate this, we examined the effect of shifting the wavelength 
grid by —0.2, —0.1, 0.1 and 0.2 pixels. The shifts this induced in 
A/i/ (j, from the orig inal value are -1.3 x 10~ 6 , -2.1 x 10" 6 , 
+ 1.0 x 10~ 6 and +0.2 x 10~ 6 . The standard deviation of these 
values is w 1.4 x 10~ 6 , and so we adopt 1.4 x 10 -6 as an estimate 
of the potential error in Ap/ p on account of the re-dispersion of 
the contributing exposures. 



6 RESULT INCLUDING SYSTEMATIC ERRORS 

In Table [3] we accumulate the potential systematic errors from our 
discussion above and give our final estimate of Ap/ p including the 
systematic component. Although the distribution of systematic er- 
rors is likely to be Gaussian in many cases, the impact on Ap,/ p 
arising from the distortion in the wavelength scale (for instance) is 
an upper limit. The probability distribution of the sum of random 



variables is given by the convolution of their individual probability 
density functions. Thus, to estimate our final uncertainty, we con- 
volve the distributions assumed for each of the sources of uncer- 
tainty, and give the standard deviation of the resultant distribution 
as our uncertainty estimate. 

This yields our final estimate of Ap/ p for the z — 2.811 
absorber toward Q0528— 250 as 

Ap 

A* 



= (0.3±3.2 stat ±1.9 s 



x 10' 



(5) 
(6) 



= (0.3±3.7)tot x 10~ 6 . 

7 OTHER ANALYSIS 

7.1 H 2 /HD column density ratio 

The HD/H2 column density ratio is of astrophysical interest. We 
calculate the total H 2 column density under the 3-parameter, b — 
F(J) model via the method given in the caption to Table [6] to be 
log 10 (iV/cm~ 2 ) = 16.556 ± 0.024. Column densities for each of 
the H2 J-levels are given in Table [6] The HD J = column den- 
sity is log 10 (A7cm- 2 ) = 13.267 ± 0.072 directly from the fit. 
We assume that the total H2/HD column density ratio is approxi- 
mated by the total H2/HD( J = 0) column density ratio ( jTumlinsori] 
|et al.|20f0) . This directly gives the HD/H2 column density ratio as 
(5.4+1.1) x 10^orlog 10 [JV(HD)/jV(H 2 )] = -3.3 + 0.2. Al- 
ternatively, log 10 [iV(HD)/2iV(H2)] = -3.6+0.2. This compares 
with log 10 [iV(HD)/2iV(H2)] = -4.1 ± 0.2 from the z = 2.059 
absorber toward J2123-0050 and -4.8 ± 1.5 from the z = 2.627 
absorber toward FJ0812+32B (Tumli nson et al.|2 010). The result 
presented here for log 10 [AT(HD) /2JV(H 2 )] from Q0528-250 dif- 
fers from the value from J2123— 0050 at the 1.8<r level, which is 
large but not grossly inconsistent, particularly since the estimate 
here is statistical only and does not consider gas cloud kinematics 
or potential systematic effects. Additionally, because the low- J H2 
transitions in the Q0528— 250 absorber are saturated it is difficult 
to accurately estimate the column density for H2 in Q0528— 250; 
|Tumlinson et al.| ( f2010| > noted this for the absorber in J2123— 0050. 
As such, some caution is warranted in comparing the HD/H2 ra- 
tio presented here to other results. Additionally, our fits were con- 
structed with the H2 line intensities fitted as free parameters. This 
makes it challenging to accurately estimate the total H2 column 
density; to obtain a more accurate HD/H2 column density ratio, 
this absorber should be re-analysed specifically for this purpose. 



8 DISCUSSION 

Firstly, it is clear that the result obtained here is consistent with 
no variation of p. Secondly, the result is consistent with the 
Q0528-250 result in |King et al] (2008) , and also the Ap/p es- 
timates obtained from Q0405-443 and Q0347-383 in |King et af] 
p008} and from J2 123 -0050 in |Malec et aT] ( |20T0t . From con- 
sideration of all the H2/HD measurements from quasar absorbers 
(which are all at z > 2) there is no statistically significant evidence 
for cosmological variation in p. In particular, we can formally com- 
bine the results from |King et al.H2008) , |Malec et al. 1(20 10) and 
this work under a weighted mean model (assuming that a weighted 
mean model is a legitimate description of the data), and we obtain 
Ap/p = (2.3 ± 2.2) x 10~ 6 , with \l about the weighted mean 



of « 0.9. The fact that ^ w 1 suggests that the uncertainty esti- 
mates on Ap/ p are approximately correct. We show the result of 
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Table 5. Error budget for the analysis of the z = 2.811 absorber toward Q0528— 250 from the sources described in section|5] The third column gives the 
magnitude of the uncertainty estimate. The final error estimate is calculated as the standard deviation of the convolution of the assumed distributions of the 
individual error estimates; the assumed distributions are given in the fourth column. 



Source of error 


A(i//i(10- 6 ) 


<5(A/Vm) (10- 6 ) 


Assumed distribution 


Systematic distortion of ThAr wavelength scale 




±1.1 


Uniform 


Intra-order wavelength scale distortions 




±0.3 


Gaussian 


Velocity structure and spatial segregation 




±0.8 


Gaussian 


Re-dispersion of spectra 




±1.4 


Gaussian 


Total systematic error 




±1.9 




Statistical 


0.3 


±3.2 


Gaussian 


Final estimate 


0.3 


±3.7 





Table 6. Column densities for the different J-levels under the 3-component, 
b = F(J) model, n gives the number of transitions contributing to a par- 
ticular J-level. Because the column densities for each transition are fitted 
as free parameters, the log 10 column density for each J-level is calculated 
as the weighted mean of the log 1( ) column densities for the transitions in 
that J-level, but where statistical errors are increased in quadrature with an 
additional term, <tj, which is calculated so that xt = 1 about the weighted 
mean. Without the inclusion of aj, xt 3> 1, which reflects the accumula- 
tion of errors in the local model for the continuum and weak, unmodelled 
Lyman-a transitions which overlap with the transitions in question, crj is 
given in the fourth column. 



J 


n 


log 10 (JV/cm- 2 ) 


lo gio(°"./) 





3 


15.68 ±0.25 


0.43 


1 


8 


16.028 ± 0.094 


0.16 


2 


24 


16.008 ±0.055 


0.23 


3 


28 


15.917 ± 0.030 


0.14 


4 


13 


14.476 ± 0.025 


0.084 


Total 


76 


16.556 ± 0.024 


n/a 



this work in comparison with other quasar constraints on A^i//x in 
Fig.© 



9 CONCLUSION 

We have used new observations of the molecular absorber at 
z — 2.811 towards Q0528— 250 to derive a new, strong constraint 
on Afi/fi at high redshift, namely Afi/n = (0.3 ± 3.2 sta t ± 
1.9 S ys) x 10~ 6 = (0.3 ±3.7) x 10" 6 . This result is the strongest 
individual constraint on variation in n at z > 1, and together with 
other measurements of H 2 /HD absorbers provides a very strong 
constraint on evolution in over most of the observable history of 
the universe. The value of Afi/fi presented here is consistent with 
the values presented in |King et aT1(2008F and |Malec et al.H2010f . 

Our result was derived using 76 H2 transitions and seven 
HD transitions. This is the first time that HD has been detected 
in this absorber, with a column density of log 10 (iV/cm - ) = 
13.27 ± 0.07. We measure the H2/HD column density ratio to be 
log 10 [jV(HD) /7V(H 2 )] = -3.3 ± 0.2 

Similar to |King et al]f2008) and |Malec et al^OlO) , we ap- 
plied the comprehensive fitting method (CFM), where we simulta- 
neously fitted the H2/HD transitions with the surrounding Lyman- 
a forest lines. By modelling the forest simultaneously with the 
H2/HD transitions, our estimate of Afi/fi should be more accu- 
rate. Similarly, the uncertainty in determining the locations of the 
forest lines directly feeds into our uncertainty estimate on A[j,/fi, 
helping to ensure that our uncertainty is not under-estimated. 

On the basis of the AICC, we concluded that a model of the 
absorber with 3 velocity components is preferred. We imposed a 
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Figure 6. Recent cosmological constraints on Afi/fi derived from quasar 
absorption measurements. Filled circles indicate the H2 measurements from 
Q0405-443, Q0347-383 and Q0528-250 in |King et al.||2008), the square 
indicates the H2/HD measurements from J2 123— 0050 in |Malec et al.| 
J2010}, the star indicates the result of this paper, the open circle indicates the 
Q0347-383 measurement of |Wendt & MolaroH2011> [which is not inde- 
pendent of the Q0347— 383 measurement in |Kmgetaf (2008 1 due to the use 
of common spectra] and the triangles indicate the constraints from ammonia 
in |Murphy etaf1(2008at , |Henkel et al.H2009) and |Kanekar|pOTT) . The two 
Q0528-250 points and the two Q0347-383 points have been slightly dis- 
placed in redshift for clarity. Where uncertainty estimates were given with 
both a statistical and systematic component, we added these in quadrature 
to derive the error bars shown. The red, solid line indicates the weighted 
mean of the z > 1 points excluding the Q0347— 383 point from Wendt & 
MolaroH201 1) [as it is not independent of the Q0347— 383 point from |King 



et al. 1(200 8'!], and the blue, dashed line indicates the la uncertainty around 
that weighted mean. Clearly, if the z < 1 measurements were included in 
the weighted mean then they would dominate the result due to their small 
error bars. 



number of assumptions about the velocity structure of the absorber 
in order to reduce the number of free parameters in the fit. Inves- 
tigation of potential differences between 'cold' (J 6 [0, 1]) and 
'warm' (J £ [2,4]) components showed that the potential impact 
of our assumptions on the estimate of Afi/fi is significantly smaller 
than the statistical error on A/j,/(i. 

We investigated a number of potential systematic effects, in- 
cluding: known wavelength calibration errors due to uncertainties 
in the ThAr calibration; observed intra-echelle order wavelength 
scale distortions of unknown origin; assumptions about the velocity 
structure, and; the effect of re-dispersing the individual exposures 
onto a common wavelength grid. The contribution of each of these 
effects is small, and in aggregate the total uncertainty in Afi/fi is 
dominated by statistical rather than systematic error sources. This 
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is reassuring, and means that future VLT/UVES observations of 
H2/HD absorbers may continue to increase the constraint onA(j//i 
in a meaninful way, provided that appropriate absorbers can be 
identified. 

Although there is no significant evidence for evolution in /x, 
and little (~ 2a) evidence for dipolar variation in fi across the sky, 
it is intriguing that the dipole direction in a dipole+monopole model 
fitted to the z > 1 H2/HD A/x//x constraints points in a similar 
direction to the Aa/a dipole described in King et al. \201 1| . Given 
the distribution of existing H2 results on the sky, new H2 absorbers 
in judicious locations will be able to place strong constraints on 
spatial variation in fj,. 

The current number of quasars known to contain molecular 
hydrogen absorption is small, and almost all existing analyses have 
focused on the objects Q0405-443, Q0347-383, Q0528-250 and 
J2123— 0050. Given that analysis of molecular hydrogen currently 
remains the most direct and precise way of investigating potential 
evolution in n at z > 2, it would be advantageous if many more 
H2/HD systems could be discovered. Although new observations 
of existing H2/HD systems can increase the SNR, common system- 
atics may remain. The use of more absorbers in different quasars 
would increase both the redshift and spatial coverage of the A^i//i 
constraints. 
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APPENDIX A: VOIGT PROFILE FITS 

Online only: Figures [AT] through |A14| show our Voigt profile model 
for the z = 2.811 absorber toward Q0528— 250 and the surround- 
ing Lyman-a forest regions, indicating both the positions of the 
H2/HD components as well as the H I components used to fit the 
surrounding Lyman-a forest. 
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Figure Al. H2/HD fit for the z = 2.811 absorber toward Q0528— 250 (part 1). The vertical axis shows normalised flux. The model fitted to the spectrum is 
shown in green. Red tick marks indicate the position of H2/HD components, whilst blue tick marks indicate the position of blending transitions (presumed to 
be Lyman-a). Normalised residuals (i.e. [data - model]/error) are plotted above the spectrum between the orange bands, which represent ±lcr. Labels for the 
H2 transitions are plotted below the data. 
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Figure A2. H2/HD fit for the z = 2.811 absorber toward Q0528— 250 (part 2). The vertical axis shows normalised flux. The model fitted to the spectrum is 
shown in green. Red tick marks indicate the position of H2/HD components, whilst blue tick marks indicate the position of blending transitions (presumed to 
be Lyman-a). Normalised residuals (i.e. [data - model]/error) are plotted above the spectrum between the orange bands, which represent ±lcr. Labels for the 
H2 transitions are plotted below the data. 
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Figure A3. H2/HD fit for the z = 2.811 absorber toward Q0528— 250 (part 3). The vertical axis shows normalised flux. The model fitted to the spectrum is 
shown in green. Red tick marks indicate the position of H2/HD components, whilst blue tick marks indicate the position of blending transitions (presumed to 
be Lyman-a). Normalised residuals (i.e. [data - model]/error) are plotted above the spectrum between the orange bands, which represent ±lcr. Labels for the 
H2 transitions are plotted below the data. 
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Figure A4. H2/HD fit for the z = 2.811 absorber toward Q0528— 250 (part 4). The vertical axis shows normalised flux. The model fitted to the spectrum is 
shown in green. Red tick marks indicate the position of H2/HD components, whilst blue tick marks indicate the position of blending transitions (presumed to 
be Lyman-a). Normalised residuals (i.e. [data - model]/error) are plotted above the spectrum between the orange bands, which represent ±lcr. Labels for the 
H2 transitions are plotted below the data. 
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Figure A5. H2/HD fit for the z = 2.811 absorber toward Q0528— 250 (part 5). The vertical axis shows normalised flux. The model fitted to the spectrum is 
shown in green. Red tick marks indicate the position of H2/HD components, whilst blue tick marks indicate the position of blending transitions (presumed to 
be Lyman-a). Normalised residuals (i.e. [data - model]/error) are plotted above the spectrum between the orange bands, which represent ±lcr. Labels for the 
H2 transitions are plotted below the data. 
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Figure A6. H2/HD fit for the z = 2.811 absorber toward Q0528— 250 (part 6). The vertical axis shows normalised flux. The model fitted to the spectrum is 
shown in green. Red tick marks indicate the position of H2/HD components, whilst blue tick marks indicate the position of blending transitions (presumed to 
be Lyman-a). Normalised residuals (i.e. [data - model]/error) are plotted above the spectrum between the orange bands, which represent ±lcr. Labels for the 
H2 transitions are plotted below the data. 
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Figure A8. H2/HD fit for the z = 2.811 absorber toward Q0528— 250 (part 8). The vertical axis shows normalised flux. The model fitted to the spectrum is 
shown in green. Red tick marks indicate the position of H2/HD components, whilst blue tick marks indicate the position of blending transitions (presumed to 
be Lyman-a). Normalised residuals (i.e. [data - model]/error) are plotted above the spectrum between the orange bands, which represent ±lcr. Labels for the 
H2 transitions are plotted below the data. 
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Figure A9. H2/HD fit for the z = 2.811 absorber toward Q0528— 250 (part 9). The vertical axis shows normalised flux. The model fitted to the spectrum is 
shown in green. Red tick marks indicate the position of H2/HD components, whilst blue tick marks indicate the position of blending transitions (presumed to 
be Lyman-a). Normalised residuals (i.e. [data - model]/error) are plotted above the spectrum between the orange bands, which represent ±lcr. Labels for the 
H2 transitions are plotted below the data. 



26 J. A. King et al. 



O - 



3980 



o - 



3990 



X 



V) 

E 



o - 



4000 



o - 




4010 



4015 



4020 



1 1 1 1 I ^J- 

_ 


II 1 1 


1 


L4P3 V . 

i 



4020 



4025 

Wavelength (A) 



4030 



Figure A10. H2/HD fit for the z = 2.811 absorber toward Q0528— 250 (part 10). The vertical axis shows normalised flux. The model fitted to the spectrum 
is shown in green. Red tick marks indicate the position of H2/HD components, whilst blue tick marks indicate the position of blending transitions (presumed 
to be Lyman-a). Normalised residuals (i.e. [data - model]/error) are plotted above the spectrum between the orange bands, which represent ±lcr. Labels for 
the H2 transitions are plotted below the data. 
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Figure All. H2/HD fit for the z = 2.811 absorber toward Q0528— 250 (part 11). The vertical axis shows normalised flux. The model fitted to the spectrum 
is shown in green. Red tick marks indicate the position of H2/HD components, whilst blue tick marks indicate the position of blending transitions (presumed 
to be Lyman-a). Normalised residuals (i.e. [data - model]/error) are plotted above the spectrum between the orange bands, which represent ±lcr. Labels for 
the H2 transitions are plotted below the data. 
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Figure A12. H2/HD fit for the z = 2.811 absorber toward Q0528— 250 (part 12). The vertical axis shows normalised flux. The model fitted to the spectrum 
is shown in green. Red tick marks indicate the position of H2/HD components, whilst blue tick marks indicate the position of blending transitions (presumed 
to be Lyman-a). Normalised residuals (i.e. [data - model]/error) are plotted above the spectrum between the orange bands, which represent ±lcr. Labels for 
the H2 transitions are plotted below the data. 



New constraint on ^-variation from Q0528—250 29 




4130 



4135 



4140 



o - 



4140 




4145 



4150 




4150 



4155 



4160 



j_ im i 



ii I ii i i 



.u 



vir 



_i i i i_ 



01 



_i i i i_ 



r- CM 



4160 



4165 



4170 



o * 




4170 



4175 

Wavelength (A) 



4180 



Figure A13. H2/HD fit for the z = 2.811 absorber toward Q0528— 250 (part 13). The vertical axis shows normalised flux. The model fitted to the spectrum 
is shown in green. Red tick marks indicate the position of H2/HD components, whilst blue tick marks indicate the position of blending transitions (presumed 
to be Lyman-a). Normalised residuals (i.e. [data - model]/error) are plotted above the spectrum between the orange bands, which represent ±1<t, Labels for 
the H2 transitions are plotted below the data. 
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Figure A14. H2/HD fit for the z = 2.811 absorber toward Q0528— 250 (part 14). The vertical axis shows normalised flux. The model fitted to the spectrum 
is shown in green. Red tick marks indicate the position of H2/HD components, whilst blue tick marks indicate the position of blending transitions (presumed 
to be Lyman-a). Normalised residuals (i.e. [data - model]/error) are plotted above the spectrum between the orange bands, which represent ±lcr. Labels for 
the H2 transitions are plotted below the data. 



